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ABSTRACT 

Initially cold and spherically symmetric self-gravitating systems may give rise to a 
virial equilibrium state which is far from spherically symmetric, and typically triaxial. 
We focus here on how the degree of symmetry breaking in the final state depends 
on the initial density profile. We note that the most asymmetric structures result 
when, during the collapse phase, there is a strong injection of energy preferentially 
into the particles which are localized initially in the outer shells. These particles are 
still collapsing when the others, initially located in the inner part, are already re¬ 
expanding; the motion of particles in a time varying potential allow them to gain 
kinetic energy — in some cases enough to be ejected from the system. We show 
that this mechanism of energy gain amplifies the initial small deviations from perfect 
spherical symmetry due to finite N fluctuations. This amplification is more efficient 
when the initial density profile depends on radius, because particles have a greater 
spread of fall times compared to a uniform density profile, for which very close to 
symmetric final states are obtained. These effects lead to a distinctive correlation of 
the orientation of the final structure with the distribution of ejected mass, and also 
with the initial (very small) angular fluctuations. 

Key words: Cosmological structure formation, gravitational clustering, N -body sim¬ 
ulation 


1 INTRODUCTION 


That self-gravitating systems initially in highly spherically 
symmetric configurations can relax to virial equilibria 
which break this sym metry strong l y has been k nown 
for several decades dPolvachenko fe Shukhmanl Il98ll: 

I Merritt fe Aguilar! 119851 ) and d ocumented since then b y 
many numerical s tu dies ( se e e.g. Aguilnj_^^Ierritt| ( 1990j); 
Theis fe Spurzeml (ll999l); IBoilv fe' Athanassoulal 1 20061) : 
Barnes. Lanzel fe Williams 1 20091) : IWorrakitpoonpoJ 
1 2014) 1. This phenomenon 
gued to play a crucial role 


This 
play a 

formation (see 


US 


of formation, 
cosmological 


and ar- 

_st ructur e 

Huss. Jain fe Steinmetd dl999l) : 


iMacMillan, Widrow fe Henriksenl ( 20061 11 has come to be 


referred to as “radial orbit instability” (ROI). This name 
has been adopted since such an inst ability has been shown 
llAntonovlll96ll : [Fridman et al.lll984 to characterize spher¬ 
ically symmetric stationary solutions of the collisionless 
Boltzmann equation with purely radial orbits. Further it i s 
plausible, as argued originally bv lMerritt fe Aguilar! (Il985ll . 
that a similar mechanism is responsible for the formation of 
triaxial structures observed starting from very cold initial 


conditions, as in this case collapse tends to produce strongly 
radial orbits. Different authors (see references above) have 
discussed how the symmetry breaking develops during the 
evolution from both si mple power law density profiles (e.g. 
IBoilv fe Athanassoula |2006|)) and from cosmolog i cal ini tial 
conditions (e.g. IMacMillan, Widrow fe Henriksen (120061) . 


In this paper we consider how the degree of the final 
symmetry breaking is related to the initial condition 
specifically to the exponent of the initial density profile 
for the case of completely cold initial conditions. Our fo¬ 
cus on this aspect of the problem allows us to elucidate the 
mechanism by which the symmetry breaking actually oc¬ 
curs in the process of collapse from cold initial conditions. 
More specifically, we show in detail how fluctuations break¬ 
ing spherical symmetry may be amplified by the very large 
energy changes characteristic of the very violent relaxation 
from cold initial conditions. This amplification is most ef¬ 
fective when the energy change a particle undergoes is both 
large and strongly correlated with its initial radial position, 
leading to a maximal effect from density profiles with in¬ 
termediate exponents. We underline that the mechanism we 
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identify bringing about the amplification of the symmetry 
breaking operates far from equilibrium and has no apparent 
link to the ROI of equilibrium systems. 


2 NUMERICAL SIMULATIONS 


For our study we have simulated numerically, using the N- 
body code Gadget dSpringelll2005l l. the evolution from initial 
conditions in which N particles are distributed randomly 
inside a sphere following a radial density profile p(r) oc r~ a , 
with a in the range 0 ^ a ^ 2.5 (the reasons for the choice 
of this range will be discussed below) The family of initial 
conditions are thus characterized by the two parameters a 
and N. We will focus here on the dependence on a of the 
degree of symmetry breaking of the relaxed state. We have 
also varied N systematically (for each a) in a range from 
a few thousand to one hundred thousand particles, and in 
this range of N our essential results and analysis are weakly 
sensitive to this parameter. We will report in future work a 
more detailed investigation of the subtle (and numerically 
challenging) issue of the asym ptotic large N dependence of 
spherical symmetry bre aking llBoilv Athanassoulal [20061 : 
IWorrakitpoonponll2014l l. 

All results presented here are for simulations in which 
energy was conserved to within a tenth of a percent. For 
simulations with a in the range [0.25, 2] this level of energy 
conservation has been attained using typical values of the 
essential numerical parameters in the GADGET code [0.025 
for the rj parameter controlling the time-step, and a force ac¬ 
curacy of of = 0.001]. The cases in which a is outside this 
range are numerically more challenging because of singular¬ 
ities — discussed further below — both for a = 0 and a = 3 
in the limit N —> oo. For these cases we have subjected 
our results to additional tests of their robustness, check¬ 
ing their stability in particular to smaller time-steps (se e 
also the discussion inljo vce. M arcos fe Svlos Labinil (l2009h : 
iBenhaiem fe Svlos Labini ( 20151 11. We have also studied 
carefully the effects of varying the force smoothing parame¬ 
ter (corresponding to its small scale at small distances), and 
we have found our results to be stable provided it is sig¬ 
nificantly smaller than the minimal characteristic size (see 
below) attained by the structure during its collapse For the 
simulations reported below the smoothing parameter is al¬ 
ways in this latter range. 


3 RESULTS 

3.1 Collapse and Virialization 


Before turning to the central issue in this article — the 
dependence of the asymmetry of the final virialized state 
on the exponent a — we consider first how various spher¬ 
ically averaged indicators of the global evolution of the 
system vary with a, and how these dependencies can 
be understood. The results in this section are in line 
with numerous previous studies of cold systems with ini¬ 
tial conditions in the class we consider, or simil a r one s 
(e.g. van Albada lil982h: Aarseth, Lin fc Papaloizoul d 1988T) : 
Ijovce, Marcos fc Svlos Labinil ( 2001^ 1. We simply focus on 
the quantities and behaviors which will be most relevant for 
our analysis of the symmetry breaking in the next section. 



Figure 1. Time evolution of the gravitational radius (as defined 
in Eq[T} normalized to its initial value R g ( 0), for different a and 
with N = 10 4 . The inset shows the characteristic time of the 
spread in fall time AT (normalized to r c ) in simulations with 
the indicated IV; the dotted line corresponds to the analytical 
prediction given by Eq. ®. 


In line with previous studies of such initial conditions, 
we observe that the system evolves through a strong collapse 
followed by a re-expansion which very rapidly leads to a 
virial equilibrium in which most of the initial mass is bound. 
FigOshows, for different indicated values of a, the temporal 
evolution of the gravitational radius defined as 


Rg(t) 


GM b (t) 

\W b (t)\ 


( 1 ) 


where M b (t) and W b (t) are respectively the mass which 
is bound (i.e. particles with negative energy) and the po¬ 
tential energy of this mass. The unit of time here is r c = 
\J 32 GPq( rIo ’ w here p 0 (Ro ) is the average mass density in¬ 
side the radius Ro of the initial spherical configuration, of 
radius Ro- It corresponds to the time for a particle initially 
at the outer periphery (i.e. at r = Ro) to fall to the center, 
in the continuum approximation (i.e. taking N —>• oo keep¬ 
ing the initial mass density profile fixed) and without shell 
crossing. The scale R g (t) is a measure of the characteristic 
size of the system, and we observe in all cases, to a first 
approximation, the same qualitative behavior — a collapse 
to a minimal size attained around t = r c followed by a re¬ 
expansion and stabilization. In all cases the bound particles 
form a virialized structure, with the virial ratio (which we 
do not display here) showing a very similar qualitative be¬ 
haviour to that of R g (t) but with a final value close to —1 
in all cases. There is, however, also a clear trend with a: 
the smaller is a, i.e. the closer to flat the density, the more 
violent is the collapse, with the system reaching a deeper 
minimum in a shorter time Q. Further, we note that the 


1 We do not show data for a = 2.5 in Fig|T]because the potential 
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larger is a the denser are the inner shells of the cloud and 
the sooner the collapse starts. 

The variation of the characteristic time for the collapse 
and re-expansion with a is quantified in the inset in the fig¬ 
ure. It shows, as a function of a, the measured time AT 
estimated as the difference between the two times at which 
R g (t) = R* g defined as R* = ( R a g syn + R™ n )/2 where R a g syn 
is the estimated asymptotic value of R g at t r c . The 
continuous curve, which has been extended to a = 3, is ob¬ 
tained from the following simple considerations. We work 
in the approximation that departure from spherical symme¬ 
try, and also the effects of shell crossing, can be neglected. 
For an initial mass density with radial profile p(r) oc r~ a 
for r < Ro, and zero for r > -Ro, mass at an initial radial 
distance ro from the center will then fall to the center in 
a time r c (r 0 ) = = T c (r 0 /R 0 ) a/2 , where p 0 (r 0 ) 

is the initial mass density of the sphere of radius ro. The 
distribution h(r) of these fall times to the origin can then 
be calculated using 47 rp(ro)rodro = Mh(r)d,T (where M is 
the total mass). One finds 


h{r) 


2(3 - a) /jr \ 3( “ _1) 
UT C \T C J 


( 2 ) 


for r ^ r c (and h(r) = 0 otherwise). The spread in the fall 
times can be characterized by the variance of h(r), 

A T th = 2V(r 2 ) - <r>2 = . (3) 

This expression, plotted in the inset of Fig. [T| reaches a 
maximum at a « 2.3, and goes to zero for both a = 0 and 
a = 3 . For a = 0, all particles fall to the origin at r c (the 
well-known singularity of the canonical “spherical collapse 
model”), while as a steepens towards a = 3 almost all of 
the mass is initially at small radii with fall times which are 
very small compared to r c . In the inset of Fig.[T| we see that 
Eq. J3| traces well the behavior of the measured AT up to 
a ~ 2. Thus, up to this value, the characteristic time of the 
variation of the total potential indeed just reflects the spread 
in the particles’ fall times. Further for larger a the behavior 
of h(r) and A Tth indeed reflect the qualitative change in 
behavior of the collapse we observe in our simulations: while 
the collapse is completed only at t ~ r c when the outermost 
mass falls, most of the mass falls at very much shorter times. 
It is for this same reason that accurate numerical integration 
becomes more costly as a increases and we report results 
only up to a = 2.5. 

One other macroscopic feature of relaxation from cold 
initial conditions of this kind, which will be relevant in our 
discussion below, is that they often lead to mass ejection 
i.e. some particles gain enough energy so that their total 
energy i s positive and they can escape to infinity. We show 
(see also lSvlos Labinil d2013h j in Fig.[5]the fraction pf of the 
particles with positive energy after relaxation (at t ~ 5r c ), 
for different a (and N). The observed behavior as a function 
of a — maximal ejection at a = 0, followed by a monotonic 
decrease (approximately exponential) with a in a range until 
a ~ 2, beyond which there is a much sharper drop — is 


energy Wb diverges for a ^ 2.5 in the limit N —► oo. As a measure 
of the characteristic size of the system we have used in this case 
the radius containing 90% of the mass. 



Figure 2. Fraction of ejected particles for different a and different 
number of particles (see labels): average and standard deviation 
over 20 realizations for N = 10 3 and N = 10 4 , and 5 realizations 
in the other two cases. 


clearly related to the qualitative behavior of the fall times 
discussed above. This will be seen more explicitly in our 
analysis below. 


3.2 Symmetry breaking of relaxed state 


We focus now on the question of how the symmetry breaking 
in the final state depends on the exponent a characterizing 
the density profile of the initial condition. Shown in the up¬ 
per plots of Fig. [3] are a projection in a chosen plane of the 
resulting virialized configurations for the various indicated 
values of a, for simulations with N = 10 4 . Visual inspection 
suggests that the breaking of spherical symmetry is appar¬ 
ently strongest in the intermediate values of a, and weakest 
for the case a = 0. This is confirmed by Fig. [3] (lower panel) 
which shows, as a function of a, the parameter iso (the “flat¬ 
tening ratio”) of the relaxed state defined as 


Ai , 
tP = A3 ~ 1 


( 4 ) 


where Ai and A 3 are, respectively, the largest and smallest 
of the three eigenvalues of the moment of inertia, and the 
subscript indicates that the estimate is made on the P % of 
particles which are most bound (we take here P = 80 fol¬ 
lowing common practice in the literature). We do not show 
here any information about the intermediate eigenvalue A 2 , 
but analysis of it shows that it is typically not close to 
the value of either of Ai or A 3 , i.e. the relaxed structures 
are quite triaxial. For each a the different points in the 
lower panel in Fig. [3] correspond to the indicated values of 
N, and the error bars to the standard deviation measured 
over the indicated number of realizations in each case. The 
plot shows clearly — in agree ment with previous studies 
of the se initia l conditions (e.g. Aarseth, Lin fe P apaloizoul 
'jovce. Marcos fc Svlos Labinil d2009h : 


iBames. Lanzel fe Williams! i 2009h : IWorrakitpoonpoa 

I 20l4) ~) — that the final state in the case a = 0 is in fact 
very close to spherically symmetric. Further we note that 
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Figure 3. Upper plots: Projections of the virialized structure for 
a realization with N = 10 4 for different a at the final time t ps 5r c . 
Lower plot: flattening ratio (average and standard deviation) tgo 
as a function of a, estimated over 20 realizations for N = 10 3 and 
N = 10 4 , and 5 realizations in the other two cases. 



Figure 4. Change in particle energies (average in bins) in units 
of eo = GN 2 m/Ro, plotted against their initial radial position ro 
(normalized to the initial maximum radius Ro), for the different 
indicated a, and N = 10 4 . Note that the scale on the y-axis is 
different in the different plots. 


above a ~ 1.5 there is a clear trend towards progressively 
weaker symmetry breaking, albeit with a lesser suppression 
than in the case a = 0. We have checked and confirmed that 
these trends of tgo with a are also observed with different 
values of P in the calculation of tp. 

In the light of the discussion in the previous section, 
these lots suggest a possible correlation between the ob¬ 
served asymmetry and the qualitative behavior of the spread 
in the fall times. Further the (relative) suppression at larger 
a suggests that there may be some connection between the 
ejection of matter and the degree of symmetry breaking. 

How are particles’ fall times related to the final state, 
and in particular its asymmetry? In the relaxation process 
from such cold initial conditions, particles energies change 
greatly as they move in the time dependent mean field. In¬ 
deed this is the essence of “viole nt relaxation” as originally 
described by |Lvnden-Belll dl967f ). It can be verified directly, 
by tracking the energy of individual particles, that the en¬ 
ergy change of any given particle occurs essentially as it 
passes through the center of the collapsing/re-expanding 
structure. This is the case simply because the mean field 
is most intense and most rapidly varying at this time. As a 
consequence this energy change depends essentially on the 
time window in which the particle passes through this cen¬ 
tral region. The correlation between initial radial position 
and fall time — which is strong except in the limit a = 0 


where all particles have the same fall time (modulo finite N 
fluctuations) — might thus be expected to lead to a corre¬ 
lation between the energy change of a particle and its initial 
radial position. 

To see whether this is indeed the case we plot in Fig|4] 
the total energy change, i.e. the difference between initial 
energy and the energy measured at a time well after the col¬ 
lapse, when the stationary state has been reached, averaged 
in bins, as a function of their initial radial position, for sim¬ 
ulations with N = 10 4 . As anticipated we see that, except 
for the case a = 0, there is indeed a very clearly identifi¬ 
able correlation between energy change and initial position: 
up to a = 2 a large positive energy boost is obtained by a 
large fraction of the particles in the outer shells, while at the 
larger values of a large energy decreases are experienced by 
the particles in the inner shells. The explanation for these 
features, and for the behavior in the case a = 0, is closely 
related to our discussion of the previous section. As noted, 
particles’ energies change essentially as they pass through 
the center of the structure, and what determines the energy 
change is the temporal variation of the mean field they move 
in at this time. 

The behavior of the mean-field potential at any point in 
the center of the structure reflect approximately that of the 
total potential shown in Fig[T| Particles which pass through 
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Figure 5. Behavior of /-soft) (circles) and /. i qq (squares, cor¬ 
responding to all the bound particles) as a function of time for 
different values of a G [0, 2.25] , in simulations with N = 10 4 . 


the center in the phase before the minimum is reached, at t ~ 
r c , will tend to lose energy because they climb out of a deeper 
potential than they fall into, while the converse will be true 
for particles which pass through the center as the system is 
re-expanding. These latter receive an energy boost, which 
has been noted to be at the origin of the mass ejection in 
both t he case of cold collapse (IJovce. Marcos fe Svlos Labinil 
120091: ISvlos Labinil I2012I. l2013h and in merging structures 
i Carucci et al J 120141) . Thus the trend we observe with a is 
clearly linked to that we discussed of the distribution of fall 
times h(r)\ for small a, a significant amount of the mass 
“falls later” and acquires an energy boost, while at larger a 
most of the mass “falls early” and loses energy. For a = 0, on 
the other hand, the correlation between the energy change 
of a particle and its initial position is very markedly weaker, 
because in this case the time of fall of a particle may become 
correlated with its initial radial position only through the 
finite N fluctuations (which regulate the singularity of the 
collapse characteristic of this case). 

Let us consider now how the large energy injection into 
the outer shells can lead to the very strong symmetry break¬ 
ing of the relaxed states observed in these cases. In the limit 
of exact spherical symmetry, the dispersion of the energy 
change at any given initial radius in Fig. Q] should vanish, 
and the observed finite dispersion is a consequence of the 
spherical symmetry breaking. Indeed what this dispersion 
implies is that the energy injection at these large radii still 
depends sensibly on the direction of arrival, and not just the 
initial radial position. 

The initial system is not perfectly spherical due to Pois¬ 
son fluctuations, and to a first approximation it can be de¬ 
scribed as an ellipsoid with i(0) ~ IV -1//2 . Therefore par¬ 
ticles experience, during the collapse, a force of which the 
tangential component depends on the angle with the ellip¬ 
soid initial major semi-axis. Consequently, as can be seen 
in Fig[5] iso (i) grows rapidly as the particles coming from 


Figure 6. Behavior of absolute value of the cosine of the angle 
between the smallest eigenvalue (which corresponds to the major 
semi axis of the ellipsoid) at time t = 0 and the smallest eigenvalue 
at time t, for different a. It can be seen that it remains close to 
unity at all times, showing that the orientation of major semi-axis 
remains essentially unchanged throughout the evolution. 


the direction orthogonal to the initial major semi-axis col¬ 
lapse to the ce nter in a shorter time tha n those along this 
axis: this is the lLin. Mestel fc Shui (Il965i ) instability whose 
effect is to amplify the initial small eccentricity. The parti¬ 
cles initially in the outermost shells which arrive latest then 
travel through the rapidly decreasing potential created by 
the other mass which is already re-expanding, and gain en¬ 
ergy. Note that, while for a = 0 the amplification and subse¬ 
quent decrease of the asymmetry to its final value occurs in 
a very short time around the collapse, for a > 0 the system 
shows a few oscillations before the complete relaxation. This 
behavior reflects that of the gravitational radius in Fig[T] In 
addition, we note that for a = 0 the amplification of the 
initial t(0) is much less marked compared to the cases with 
a > 0, in which the spread in fall times of the particles is 
much greater. For a ^ 2, most of mass now collapses at very 
short times and accordingly tso (t) shows a very fast growth 
for t < t c followed by a decrease when the external particles 
pass through the center. 

The hypothesis that the amplification of the initial tri- 
ax iality is due to finite N eff ects via the instability described 
bv lLin, Mestel fo Shul (1 19651 ). and the subsequent role of en¬ 
ergy gain in amplifying the triaxiality, can be tested for 
straightforwardly. Firstly, we would expect that the major 
semi-axis at any time t should be correlated to its value 
at t = 0. Secondly, when there is a significant mass ejec¬ 
tion, we would expect there to be a correlation between the 
angular distribution of the ejected particles and the orienta¬ 
tion of the final triaxial structure, and more specifically be¬ 
tween the preferred direction for ejection and the elongated 
axis of the final structure. As a measure of these spatial 
distributions we use simply the inertia tensor, determining 
its eigenvalues and eigenvectors. Fig[6] (upper panel) shows 
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Figure 7. Histogram of the values measured, in 20 realizations 
of the case a = 0.5 and N = 10 4 , of the cosine of the angle 6 
between the longest axes (determined from inertia tensor) of (i) 
the relaxed and ejected mass distributions (left panel), and (ii) 
the initial and ejected mass distributions (right panel). 


that the orientation of the major semi-axis indeed remains 
almost the same throughout the collapse and virialization. 
The left panel of Fig. [7] shows, for the case a = 0.5, which 
has both a strongly triaxial final state and significant mass 
ejection, a histogram of realizations of the modulus of the 
cosine of the angle between the eigenvectors corresponding 
to the smallest eigenvalues for the ejected mass and the 80% 
most bound mass. There is again a very clear positive signal 
for the correlation of the orientation of the axes. 

The right panel of Fig. [7] shows the correlation, mea¬ 
sured in the same way, for a number of realizations, be¬ 
tween the initial distribution and the ejected mass. We ob¬ 
serve again a very clear correlation, and a similar result 
is found considering the initial and final bound mass. The 
reason for this strong correlation of the orientation of the 
longest axes is simple: the moment of inertia of the initial 
mass gives a measure of the small effective anisotropic due 
to the finite N fluctuations, and more specifically the axis 
of the smallest eigenvalue is the axis along which the mass 
is “stretched” furthest away from the plane orthogonal to 
it passing through the center. It is precisely along such a 
“stretched” direction that one expects the particles to ar¬ 
rive slightly later than the others - just as the mass along 
the longest axes of an ellipsoidal distribution — and thus to 
receive a slightly larger energy kick. This provides further 
convincing evidence that the elongation of the final struc¬ 
ture indeed has its origin in the (relative) delay in particles’ 
fall along these directions. 

This mechanism of symmetry breaking is more efficient 
for a in the range [0.5, 1.5] (see Fig. [4j. There is indeed 
greatest symmetry breaking in the final state for interme¬ 
diate values of a where the energetically boosted particles 
are well localized in the outer radii, and where these same 
particles represent a significant fraction of the mass. The 




Figure 8. Density (upper panel) and radial velocity (bottom 
panel) profiles for the final configuration at t = 5r c and for dif¬ 
ferent values of a and with N = 10 4 . The dashed line has a r ~ 4 
behavior in the case of the density profile (upper panel) and a 
r~ 1 behavior in the case of radial velocity profile. 


strong suppression of asymmetry in the case a = 0 is, con¬ 
versely, due to the fact that, although there are many par¬ 
ticles which pick up large energy boosts, they come from 
many different parts of the initial structure. This leads to 
an effective averaging over the angular fluctuations and a 
much more spherical structure. 

We note that for a ^ 2 there is a markedly different 
behavior of iioo and iso- This is a result of the fact, as seen 
in Fig[2] that there is no mass ejection in these cases: late- 
arriving particles do not gain enough kinetic energy to es¬ 
cape from the system. In this case there are however a sig¬ 
nificant fraction of high energy (but bound) particles which 
reach large distances, giving rise to a configuration that is 
more asymmetric than that of the 80% most bound particles. 
Despite these differences the density and radial velocity pro¬ 
files are very similar for all values of a (see Fig(8j, showing 
res pectively a d e cay n (r) ~ r~ 4 and ( v%) ~ r _1 as noticed 
by lyi os Labini (l2013h . However the a = 0 case leads to a 
more compact configuration than a > 0 as a consequence of 
the larger spread of the energy distribution. 

We remark finall y that a very similar dynamical pro- 
cess was observed by iBenhaiem fc Svlos Labinil (|2015l i for 
the case of initial con ditions given by uniform ellipsoids. In 
addition, we note that lTheis fc Spurzeml dl999lf considered a 
Plummer density profile with a very small initial virial ratio, 
finding that a fast dynamical collapse generates at its end 
the maximal triaxiality of the system: the dynamical mech¬ 
anism that generates such a triaxial structure is the same at 
work in the power-law density profile. 
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4 DISCUSSION AND CONCLUSIONS 

In summary, the mechanism for the generation of the very 
strong spherical symmetry breaking observed for certain 
density profiles during violent relaxation from cold spherical 
initial conditions is essentially the existence of a preferential 
axis for the large “energy injection” to particles in the outer 
parts of the initial structure, which leads to an elongation of 
the final structure along the same axis. This axis is itself de¬ 
fined by the finite N fluctuations breaking spherical symme¬ 
try in the initial conditions, corresponding to an axis along 
which the matter is on average further from the origin. Par¬ 
ticles initially at larger radii along this axis fall through the 
collapsing region later than particles along the other axes, 
and as a consequence pick up a larger energy injection. After 
this time particle energies change negligibly, and the strong 
variation with energy as a function of angle along the (pre¬ 
dominantly radial) orbits leads, after further phase mixing, 
to a virialized structure with a spatial structure reflecting 
the energy injection. In particular the major semi-axis of 
the final structure is correlated both with the (very slightly) 
long axis of the initial condition, and with that along which 
the energy kick obtained by particles during violent relax¬ 
ation is greatest. This mechanism has no apparent relation 
to the instability of equilibrium systems with radial orbits. 
More specifically, the energy injection which is its essential 
ingredient occurs in a very short time during violent relax¬ 
ation when the system is very far from equilibrium. The sys¬ 
tem never approaches close to an equilibrium configuration 
which is spherically symmetric with purely radial orbits. 

The analysis presented here is of completely cold ini¬ 
tial conditions only. We can anticipate that both how the 
asymmetry of the final state, and the qualitative features of 
the collapse leading to it, will show the same behavior as a 
function of a for simple distributions of non-zero initial ve¬ 
locities, provided the initial viri al ratio b = 2 K/\W \ associ¬ 
ated is sufficiently small (see e.g. lSvlos Labinil J2012ll . which 
studies in particular the e nergy injection and mass ejection 
as a function of b). Giv en llPolvachenko fe Shukhmanlll98ll : 
iMerritt &: Aguilar! 1985ll that warm (and, in particular, equi¬ 
librium) initial conditions are observed in many cases to give 
rise to triaxial structures, by a mechanism which clearly is 
intimately related to the ROI of equilibrium systems, we are 
led to the conclusion that there are (at least) two distinct 
mechanisms subsumed under what is usually called ROI. We 
note that this conclusion is not only consistent with previous 
studies, but gives an explanation of one of their striking (and 
puzzling) results, namely that when the initial virial ratio b 
is varied, there is a critical value at which there is a qualita¬ 
tive change in behavior: above this value symmetry break¬ 
ing occurs only if the velocity distribution is sufficiently ra¬ 
dial, while below this value symmetry breaking occurs irre¬ 
spective of whether the velocity distribution is isotropic o r 
not (see e.g. Figure 4 of lBarnes, Lanzel fo Williams (|2009h . 
which locates the value at b m 0.05 — 0.15 depending on the 
profile). 

While the dependence on the velocity anisotropic is a 
“smoking gun” for “real” ROI operating for the warmer ini¬ 
tial conditions, the presence of a threshold below which the 
details of the velocity distribution has no relevance is ex¬ 
plained very naturally as the dominance in this region of 
a distinct mechanism operating far from equilibrium, dur¬ 


ing the coll ap se, as we have de scribed here. A recent study 
(IPakter, Marcos fe Levinll2013h performs an analysis of the 
stability to elliptical perturbations of a uniform sphere (i.e. 
a = 0 here) with an isotropic velocity distribution and oscil¬ 
lating under its mean field, and finds that there is a critical 
value of the initial virial ratio below which there is such 
an instability. Whether or not this is in fact the essential 
instability leading to deviation from spherical symmetry at 
early t imes during the collapses we ha ve studied, the anal¬ 
ysis of iPakter, Marcos fe Levfnl (l2013l l illustrates that, far 
from equilibrium, there are indeed such instabilities even 
when the velocity distribution is isotropic, and which are 
thus physically distinct from the ROI mechanism for equi¬ 
librium systems. 

We note finally that our analysis here has been for iso¬ 
lated systems with initial power law initial conditions in a 
non-expanding universe. A very similar phenomenology of 
symmetry breaking starting from cosmological type initial 
co nditions i n an expandi n g background ha s been described 
in lMacMillan, Widrow fe Henriksenl d2006l), and the RO I in 

this c ase has been linked fsee also iHuss. Jain fe Steinmetd 

(ll999lB to the generation of a “universal” NFW-type den¬ 
sity profile in this context. We do not observe the ROI in our 
study to be associated with such a final density profile: our 
profiles are well characterised by a quite flat inner core and 
power law decay ~ r~ 4 . To determine whether the specific 
mechanism of amplification of symmetry breaking we have 
described here — associated with the energy injection to 
material initially in the outer shells — is at play, or not, in 
the formation of triaxial dark matter halos in cosmological 
models would require further extensive study, incorporat¬ 
ing a careful comparison between the dynamics of isolated 
structures and halos in the cosmological setting. 

Numerical simulations have been run on the Cineca 
PLX cluster (project ISCRA QSS-SSG). In addition, this 
work was granted access to the HPC resources of The In¬ 
stitute for scientific Computing and Simulation financed by 
Region lie de France and the project Equip@Meso (reference 
ANR-10-EQPX- 29-01) overseen by the French National Re¬ 
search Agency (ANR) as part of the Investissements dAvenir 
program. 
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